Reflection and Exclusion of Shear Zones in Inhomogeneous Granular Materials 
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Shear localization in granular materials is studied experimentally and numerically. The system 
consists of two material layers with different effective frictions. The presence of the material in- 
terface leads to a special type of "total internal reflection" of the shear zone. In a wide range of 
configurations the reflection is characterized by a fixed angle which is analogous to the critical angle 
of refraction in optics. The zone leaves and reenters the high friction region at this critical angle 
and in between it stays near the interface in the low friction region. The formalism describing the 
geometry of the shear zones and that of refracted and reflected light beams is very similar. For the 
internal visualization of shear localization two independent experimental techniques were used (i) 
excavation and (ii) Magnetic Resonance Imaging. 
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I. INTRODUCTION AND MODEL 



Shearing of complex materials e.g. non-newtonian flu- 
ids, colloids, emulsions, foams or granular materials of- 
ten leads to shear localization [IHB] . The region with the 
highest shear rate is called the shear zone (or shear band) 
[7M21j. Experiments and numerical simulations for gran- 
ular materials showed that the shear zone can exhibit 
nontrivial curved or refracted shapes even for the sim- 
plest stationary case, where the position and shape of the 
shear zone remains more or less the same during the shear 
deformation [5lU3|. The case of inhomogeneous materi- 
als is especially important for understanding various nat- 
ural phenomena or industrial problems. In such systems, 
where the internal friction of the material changes with 
location, the shear is preferentially localized to regions 
with lower friction. A simple example is the case of a 
layered material where different layers correspond to dif- 
ferent friction — a configuration which is often seen in 
nature. The aim of the current paper is to characterize 
the puzzling shapes of shear zones in such systems by 
experimental and numerical techniques. 

The shape of the shear zone can be captured by a sim- 
ple variational model [22] where the basic idea is the 
following: the shear zone is modeled as an infinitely thin 
sliding surface. For any potential sliding surface that is 
consistent with the boundary conditions, one has to de- 
termine its total shear resistance against the driving. The 
actual shear zone will correspond to the sliding surface 
that has the lowest total resistance. This is the weak- 
est surface that yields first and separates the system into 
two solid blocks that slide next to each other. The proto- 
col for the determination of the shear resistance depends 



on the actual configuration. On one hand, for Couette 
or cylindrical split-bottom shear cells the potential slid- 
ing surfaces are loaded by the same mechanical torque. 
Therefore the weakest surface has to be determined based 
on the ability to transmit the torque. On the other hand, 
if straight shear cells are used as in the present study, the 
potential sliding surfaces experience an equal shear force. 
Therefore the shear zone will correspond to the potential 
sliding surface that is the weakest regarding its ability to 
transmit the total shear force between its two sides. In 
both cases - described by minimum torque or minimum 
force - the shear deformation corresponds to the lowest 
possible energy dissipation |10l E2] ■ 

In this paper we deal with straight shear cells where 
shearing is provided by the relative motion of two long 
parallel sliders (Figs. [I] and [2]). Then the shear zone is 
created between two material blocks, each stuck to one 
of the sliders. The location of the emerging shear zone 
can be determined by the above consideration. As the 
geometry is translation invariant in the direction of shear 
(x direction) the shear zone can be represented by a single 
path in the cross section of the cell (yz plane) . The total 
shear force that can be transmitted by such a path is 
proportional to the integral 



ds 



(1) 



where r, the maximum shear stress that the material 
can exert locally, is integrated along the length s of the 
path. The shear zone corresponds to the path for which 
the integral takes its minimum value. It is assumed that t 
is given by the product of the effective friction coefficient 
H and the local pressure [23] P 
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T = up. 



(2) 




FIG. 1: A gedankenexperiment: a material consisting of two layers is sheared by moving the bottom part of the cell according 
to the red arrow. The four panels illustrate how the boundary conditions affect the position of the shear zone in this closed 
system. From (a) to (d) the left boundary of the zone (position of the slider edges) is moved downwards step by step, whereas 
its right boundary is kept fixed. The case (a) is analogous to the refraction of light beams, (b) Limit angle of refraction /3 C . 
Panel (c) shows the "total reflection" of the shear zone. It is quite different from the total reflection of geometric optics as part 
of the zone is excluded from the high friction material and the zone arrives at the interface at the limit angle f3 c . (d) Straight 
zone, not influenced by the interface. 



In an inhomogeneous system the effective friction // 
can vary from place to place, e.g. when layers of different 
materials are present. When p is constant, the condition 
|l| reduces to J /i ds = min., which is formally identical 
to the well known Fermat principle in optics. According 
to this, the light beam traveling through an optically 
inhomogeneous material finds the optimal way, where the 
length weighted by the refraction index is extremal. This 
results in refraction of light beams at the boundaries of 
regions with different optical indices. The analogy leads 
to the idea that similar refraction effects are expected 
in these two distant fields of physics. The case of zone 
refraction can then be described by a law which is similar 
to the refraction law in optics 



as close as possible to the interface but in the low friction 
region. Thus this part of the zone enters with the coeffi- 
cient Hi into the calculation of the total transmitted shear 
force. The shear zone from both sides reaches the inter- 
face at the critical angle j3 c defined by sin/3 c = fii/fih- 
Interestingly this angle characterizes a wide range of con- 
figurations not only the limit case of refraction (Fig. . 

It can be easily verified that the above exclusion ef- 
fect and the piecewise straight shape of the shear zone 
(Fig. [IJ:) indeed correspond to the least possible shear 
force within the framework of the variational model of 
the infinitely narrow zone. In the following we test this 
picture experimentally. 



sin a 



Sill/ 



(3) 



where a and /3 are the two angles of incidence (see 
Fig. [Ik), while fn and /x/, stand for the effective friction 
in the two layers with low friction and high friction, re- 
spectively. This law was recently tested numerically and 
experimentally [TTJMl3| . 

In the present work we do not deal with the case of 
refraction but focus on configurations in which the zone 
after visiting the interface eventually returns to the high 
friction part of the sample. The boundary conditions, 
where both ends are located in the high friction mate- 
rial, correspond to the total internal reflection in geo- 
metric optics. However, in this situation the behavior of 
the shear zone differs from the behavior of light beams 
(Fig. [Ik). Our version of Format's principle, i.e. the 
global minimum of J /ids, suggests the following inter- 
esting scenario. Here the middle part of the shear zone is 
excluded from the high friction material and it is located 



II. EXPERIMENTAL METHODS 

For the experimental realization of the above described 
systems we use two materials with different frictional 
properties. One material is corundum which consists of 
angular grains (see Fig. |2p). Therefore it has higher ef- 
fective internal friction than the other material consist- 
ing of glass beads (Fig. [2]i). To characterize the dif- 
ference in the effective frictions, the angles of repose 
(6 r ) were determined by the method used in [23] and 
were df a = 21.9° for glass beads while 9 c r or = 33.2° for 
corundum. The latter is somewhat higher than the typ- 
ical value for sand (6 s r an = 30.5°). The effective fric- 
tions are simply estimated by fj, rs tanf? r . The ratio 
tan^°7tan^ Q = 1.63 gives a reasonable contrast for 
the effective friction. 

Two experimental geometries were investigated. In the 
first one — which is referred to as "closed system" (see 
Fig. [2k) — the zone was forced to start and end in the high 
friction part of the sample. Here the granular materials 
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were layered horizontally with the low friction material 
placed on top of the high friction material and shearing 
was obtained by a slow translation of the bottom plate 




FIG. 2: (a)-(b) Schematic illustration of the two experimental 
configurations used. The layered granular material is sheared 
by moving (a) the lower boundary for the "closed system" or 
(b) one of the L shaped walls for the "open system" according 
to the red arrow. The location of the shear zone is indicated 
with a white line. The upper surface of the "closed system" 
is loaded with a weight 3 times heavier than the granular 
material below, (c)-(d) Microscopic images of the materials 
used in the experiments (c) corundum (grain size d = 0.33 ± 
0.02 mm) and (d) glass beads (grain size d — 0.48 ±0.02 mm). 

according to the red arrow. If the thickness of the lower 
(high friction) layer is not too large the zone escapes it 
and then comes back to it, as it is indicated by the white 
line. In order to reduce the relative pressure difference 
between the top and bottom of the system the upper 
surface was loaded with a weight which was about 3 times 
heavier than the granular material below. The length of 
the cell was 70 cm and it had an internal cross section of 
4.5 cm x 4 cm. 

The second experimental apparatus referred to as 
"open system" (see Fig. |2]b) included two 70 cm long L 
shaped sliders (cell wall), one of which was slowly trans- 
lated in the experiments. When the two granular layers 
(with low and high friction) are arranged vertically (as 
seen on the the sketch), the shear zone takes a form as 
it is indicated by the white line. It is forced to start in 
the high friction part of the sample but is free to select 
the position of the other end at the top surface. This cell 
had an internal cross section of 4 cm x 3.5 cm. The to- 
tal shear displacement achieved in each experiment was 
between 5 and 6 cm for both cases (Figs. [2^ and(2j}). 

Two methods were used for the measurement of the 
deformation: (i) using colored samples and excavating 
the material layer by layer (see Fig. [3]) the displacement 
was detected optically by digital imaging and (ii) the 



whole cell was placed into a Magnetic Resonance Imaging 
(MRI) apparatus and the displacement was detected by 
tracer particles (see Fig. [I]) . In the following we demon- 
strate both methods. 

In the case of optical detection we used two samples 
with different colors each for both materials. In the upper 
layer, brown and gray correspond to glass beads (grain 




FIG. 3: Demonstration of the excavation method for the 
case of the closed system, (a) Sample images taken during 
excavation. The displacement of the grains inside the sam- 
ple is visualized by carefully removing the material layer by 
layer, (b) Displacement D(y, z). The position of the interface 
between the two materials was at height 1.0 cm. 

sizes d = 0.56 ± 0.02 mm and d = 0.48 ± 0.02 mm), 
while in the bottom layer, dark blue and white corre- 
spond to corundum (grain sizes d = 0.33 ± 0.02 mm and 
d = 0.23 ± 0.02 mm). After each experiment the ma- 
terial was carefully removed layer by layer [51 [T2] and 
the displacement profile D(y) was determined for each 
horizontal layer as it is illustrated by sample images in 
Fig. [3^,. We used a commercial vacuum cleaner with ad- 
ditional extension tubes to provide a slow flow rate. The 
procedure takes several hours for one experiment. Typi- 
cally 20-25 layers were recorded with the smallest inter- 
slice distance of 1 mm. The D(y,z) displacement profile 
shown in Fig.[3j3 resulted after combining the information 
obtained for all slices. 

Visualization of the internal deformation of the mate- 
rial by Magnetic Resonance Imaging was realized using 
a Bruker BioSpec 47/20 MRI scanner operating at 200 
MHz proton resonance frequency (4.7 T) at the Leib- 
niz Institute for Neurobiology, Magdeburg. The cross 
section of our experimental apparatus was optimized for 
the geometry of this device with 7 cm internal coil diam- 
eter. Since the granular samples used (corundum and 
glass beads) do not give an MRI signal, tracer parti- 
cles were needed. The tracer particles should be large 
enough allowing the detection of single particles, but also 
not too large for the best spatial resolution. The best 
signal was obtained using poppy seeds with the size of 
d = 0.75 ± 0.04 mm). Similarly to the case of excava- 
tion horizontal slices were obtained with the intcrslicc 
distance of 0.8 mm and the in plane resolution of 0.156 
mm/pixel. The slider was displaced by 2.5 mm between 
subsequent MRI scans, Fig. [4]: shows 3 subsequent images 
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FIG. 4: Demonstration of the measurements using Magnetic 
Resonance Imaging for the case of the open system, (a)-(b) 
The experimental cell filled with the material doped with 
poppy seeds, (c) Overlapped horizontal sample images taken 
at z — 3.2 cm. Poppy seeds appear as distinct spots on the 
binarized images with different colors for the 3 subsequent 
images. As it is expected only the right hand side of the ma- 
terial is moving (see Fig.JSJa), and the shear zone is shifted to 
the right with respect to the cell middle. 



overlapped taken at z = 3.2 cm. During one experiment 
20-30 displacement steps were recorded yielding a total 
displacement of 5 — 7.5 cm. The total measurement took 
about 4 hours. We note, that a potential segregation of 
the tracers during shear if present, would be immediately 
evident in the MRI images. Such segregation played no 
role in the present experiments. 



III. SIMULATION METHOD 

The computer simulations presented here are based 
on the variational narrow band model that we summa- 
rized in the introduction. However, the simulation tries 
to overcome one serious drawback of the narrow band 
model, namely, that it regards the shear zone infinitely 
thin. In reality the zone has a non-negligible finite width 
and the system has a smooth shear profile. In the simula- 
tion, wide shear zones are achieved by introducing fluctu- 
ations due to disorder of the granular material as follows. 
For a given state of the system, i.e. for one random re- 
alization of the material, we identify the weakest sliding 
surface according to the variational narrow band model 
(details are given below) . This procedure can be repeated 
many times for various random realizations which results 
in a large number of narrow shear bands. Taking an en- 
semble average over the random realizations provide us 
with a smooth shear profile and a wide shear zone. A 
justification of this approach can be seen in the three- 



dimensional structure of the actual shear zone. It adopts 
an averaged profile in the yz cell cross section by averag- 
ing over all slices along the x direction. 

Such a protocol has been first used in computer simu- 
lations performed by Torok et. al. [8], however, there are 
two major differences between those simulations and the 
present method. First, here we use independent random 
realizations as opposed to the self-organized random po- 
tential used earlier [S]. Second, the square lattice that 
was used in the previous studies [S| inevitably introduces 
preferred orientations for the shear zone. It is important 
that we avoid such a bias. Therefore here an isotropic 
random mesh is used instead of a regular lattice. 

The two-dimensional random mesh is generated in the 
yz cross section of the system perpendicular to the shear 




FIG. 5: The random mesh used by the MC simulation and 
one shear band. This situation corresponds to the straight 
split bottom cell filled with one material. The band starts at 
the split line of the bottom and ends at the free top surface 
of the system. The shear band represents the weakest sliding 
surface for the given random realization of the material. 

direction. First, we start with the grid points of a square 
lattice (without bonds) of lattice constant a. Second, the 
position of each point is regenerated in a square of size 
2a times 2a which is centered on the original position 
of the point. The new random position is then chosen 
uniformly inside the square. In the third step the points 
are connected with bonds using Delaunay triangulation 
(Fig. [5]). The resulting mesh is macroscopically isotropic 
to a good approximation. For example testing the bond 
orientations shows that all directions are equally probable 
and we found no sign of the orientation of the original 
square lattice. 

The triangles of the mesh correspond to microscopic 
blocks of the material while bonds represent interfaces 
between the blocks where sliding can take place. The 
maximum shear force |25j S that a given bond can exert 
without sliding is given by l^ipR, where I is the length 
of the bond, /i is the local coefficient of the effective fric- 
tion, p is the local pressure and R is a random number 
chosen uniformly between and 1. Each bond has its 
own random number R that is regenerated for each ran- 



dom realization of the material. It is assumed that the 
local pressure p is hydrostatic, i.e. p is proportional to 
the weight of the material above. 

Any potential sliding surface which is consistent with 
the boundary conditions can now be represented by a 
continuous chain of bonds. For example the two ends 
of the chain are fixed for the closed system according to 
Fig. [2^,, while for the open system (Fig.J^)) only the lower 
end is fixed at the split line of the boundary, the upper 
end is free, i.e. it can be anywhere at the top surface of 
the material. For such a chain V the total shear force 
that it can resist is given by the sum of bond forces Si 
for all the bonds i contained in T. The actual shear band 
r, i.e. the weakest sliding surface, is determined by the 
condition 



^ hpiPiRi = min. 



(4) 



The shape of the shear band depends on the actual ran- 
dom realization. Such a shear band represents a displace- 
ment jump and therefore a noncontinuous displacement 
field. It divides the material into two regions: The region 
on one side of the shear band has a unit shear displace- 
ment in x-direction while the region on the other side has 
zero displacement. 

In our Monte Carlo simulation we generate shear bands 
for many random realizations (typically a couple of thou- 
sand bands are collected) . The displacement field that is 
obtained by averaging over random realizations is contin- 
uous and can be directly compared to the experimental 
displacement profiles. 

The local friction coefficients \ii are set according to the 
experimentally measured values. There is only one free 
parameter left in the simulation, namely, the resolution 
of the mesh based on the lattice constant a. The parame- 
ter a has to be calibrated once for each type of material. 
This has been done with the help of the straight split 
bottom cell using only a single material and by match- 
ing the width of the shear zone between simulation and 
experiment. This provided us with calibration ratios be- 
tween grain size d and lattice constant a. We obtained 
a/d = 0.42 for glass beads, and a/d = 0.11 for corun- 
dum. This is in accordance with other measurements on 
homogeneous materials [5] showing that the thickness of 
the zone is larger for spherical beads than for irregular 
grains. This also means that in the simulations of the 
inhomogeneous systems presented in this paper a is not 
constant in space. The applied resolution of the under- 
lying mesh is different for the different materials. After 
this calibration is done there is no fit parameter left in 
the simulation method. 



IV. RESULTS 

We first present the results obtained for the closed sys- 
tem using excavation. The final displacement D(y, z) of 
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FIG. 6: (a)-(d) Displacement D(y, z) of the material in the 
closed system for four values of the position of the interface 
h — 0.5 cm, h = 0.8 cm, h — 1.0 cm and h = 1.4 cm mea- 
sured by excavation, (e)-(h) Results of numerical simulations 
obtained by the fluctuating narrow band model for the same 
configurations. 



the material is presented in Figs. [6^-d for four measure- 
ments in such a way, that regions with highest shear de- 
formation are visualized. As it is seen both — the experi- 
mental data and the numerical calculations (Figs.[6^-h) — 
confirm, that for a thin enough corundum layer the shear 
zone escapes the corundum and stays excluded from the 
corundum just above the interface and then returns to 
the corundum to reach the other corner of the setup. For 
a critical thickness h c of the corundum layer the expres- 
sion ([I]) has two minima of equal shear force, meaning 
that there are two optimal configurations: (i) the zone 
stays in the lower layer or (ii) the zone escapes the lower 
layer as we have seen for thin layers. These two configu- 
rations correspond to the same total force (Eq. |l])). For 
homogeneous pressure this yields 



h r = 




(5) 



where W stands for the width of the cell. For the 
given cell thickness the estimation for h c would be 1.1 
cm. A similar, but somewhat more complicated calcu- 
lation that includes hydrostatic pressure and the den- 
sity difference between corundum and glass beads yields 
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h c « 1.3 cm. For thick corundum layers the zone stays in 
the corundum, for which an example is shown in Fig. [6ji 
at h = 1.4 cm. The observations show very nice quan- 
titative agreement with the above calculations, since we 
find that when h is near the calculated critical thick- 
ness, the zone splits up into two branches as it is seen 
in Fig. [7j This happens at around h = 1.1 cm in the 
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FIG. 7: Splitting of the zone, (a) Displacement D(y,z) and 

(b) strain in the material in the closed system for h = 1.1 
cm of the position of the interface, measured by excavation. 

(c) -(d) Results of numerical simulations obtained by the fluc- 
tuating narrow band model. In the simulation the splitting 
of the shear zone occurs at a somewhat higher position of the 
interface (h — 1.25 cm). 



experiments and h = 1.25 cm in the simulations. The 
zone is remarkably bent in all cases, which is apparently 
a direct consequence of the vertical pressure gradient in- 
side the granular material. Namely, frictional forces are 
smaller at higher levels. Thus Eq. ([!]) results in a curved 
trajectory even in a homogeneous layer. 

Numerous experiments have been carried out using 
the closed system. In Fig. [8] we plot the sine of the 
angle of incidence as a function of the vertical posi- 
tion h of the interface between the two layers. As it 
is seen, the measured data match the expected value 
l/sin/3 c = fih/fk — 1-63 within experimental error. 

Experiments for the open system were carried out using 
both detection methods: excavation and MRI. Several 
configurations were prepared with different position of 
the vertically aligned interface. The results of three mea- 
surements are shown in Fig. [9] for y — 0.3 cm, y = 0.66 
cm and y = 0.94 cm. As it is seen the zone starts from 
the middle of the cell and moves up toward the right hand 
side of the sample. Then it leaves the high friction part 
and continues toward the top surface in the low friction 
part of the sample. We can calculate the critical dis- 
tance y c at which zone splitting is expected in a similar 
way as we did for the closed system. Here we also take 
into account that pressure is zero at the top surface and 
it increases linearly with depth. This yields the critical 
distance 
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FIG. 8: The sine of the angle of incidence f5 c as a function 
of the vertical position h of the interface between the two 
layers. Each measurement provides two data points, which 
are marked with different symbols for clarity. The horizon- 
tal dashed line denotes the ratio of the effective friction of 
corundum and glass beads estimated by angle of repose mea- 
surements. 



Vc = H 




(6) 



which gives y c — 0.93 for the filling height of H = 3.4 
cm. The simulation nicely reproduces the splitting of the 
zone for this configuration (see Fig. [9j) and it is partly 
visible on the experimental data (Fig. [9]:). In case the 
interface is further away from the center (y 3> y c ) we 
obtain a vertical shear zone in the middle of the cell (not 
shown here). Then the system behaves as if only the 
high friction material was present and the interface has 
no effect on the properties of the shear zone. 

Finally, we discuss one difference in the MRI and exca- 
vation techniques. While the MRI analysis provides the 
differential displacement in individual steps, the excava- 
tion yields only the integral displacement at the end of 
the experiment. If the zone position and shape is sta- 
tionary during the whole experiment, these two infor- 
mations are equivalent. However, if an initial transient 
corresponding to the shear zone formation is present, it 
affects only the excavation data. As the comparison of 
Figs. |9^-c with Figs. [9ji-f shows, this influence plays no 
role for the present experiments. Nevertheless the quan- 
titative characterization of the initial transient behavior, 
which is in principle possible with the MRI technique, 
will be an interesting aspect of future research. 



V. SUMMARY 

Shear localization in inhomogeneous granular materi- 
als has been studied experimentally and numerically. The 
shear zone preferentially develops in regions with lower 
friction. In systems consisting of layers with different ef- 
fective frictions the zone often changes direction at the 
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FIG. 9: Displacement D(y,z) of the material in the open system for three values of the position of the interface y = 0.3 cm, 
y = 0.66 cm and y — 0.94 cm measured by (a)-(c) excavation and (d)-(f) Magnetic Resonance Imaging, (g)-(i) Results of 
numerical simulations obtained by the fluctuating narrow band model for the same configurations. 



layer boundary — a phenomenon analogous to light re- 
fraction in geometrical optics. The formalism describing 
the geometry of the shear zones and that of refracted and 
reflected light beams is very similar. 

Here we have shown that total internal reflection ex- 
ists also for shear zones. However, unlike in optics the 
zone reflection occurs always at the critical angle of re- 
fraction (3 C . In case of shear zones this angle is defined by 
the ratio of the effective frictions of the two material lay- 
ers, sin/3 c = fiif/ih. This special reflection also involves 
that a part of the shear zone is trapped at the interface 
of the layers. We expect that such effects are present 
in naturally layered systems consisting of not only two 
but multiple layers, and could be best observed if the 



direction of shear velocity lies parallel while its gradi- 
ent is perpendicular to the layers. Our measurements 
also demonstrate, that Magnetic Resonance Imaging is 
a very powerful tool to accurately monitor the internal 
reorganization of granular materials in quasistatic flows. 
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